A priori Wannier functions from modified Hartree-Fock and Kohn-Sham equations 
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The Hartree-Fock equations are modified to directly yield Wannier functions following a proposal 
of Shukla et al. [Cliem. Phys. Lett. 262, 213-218 (1996)]. This approach circumvents the 
a posteriori application of the Wannier transformation to Bloch functions. I give a novel and 
rigorous derivation of the relevant equations by introducing an orthogonalizing potential to ensure 
the orthogonality among the resulting functions. The properties of these, so-called a priori Wannier 
functions, are analyzed and the relation of the modified Hartree-Fock equations to the conventional, 
Bloch-function-based equations is elucidated. It is pointed out that the modified equations offer a 
different route to maximally localized Wannier functions. Their computational solution is found to 
involve an effort that is comparable to the effort for the solution of the conventional equations. Above 
all, I show how a priori Wannier functions can be obtained by a modification of the Kohn-Sham 
equations of density-functional theory. 
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I. INTRODUCTION 

Ever since the introduction of Wannier functions 
in 1937ji they have been used as an alternative repre- 
sentation to Bloch functions in the study of crystalline 
solids They provide a local, atomic-orbital- 
like view on one-particle states of crystals, which is the 
direct generalization to periodic systems of the concept 
of localized molecular orbitals used by chemists to study 
bonding in molecules and clusters. In such a way, they 
come much closer to the intuitively accessible concepts of 
chemical bonding than the conventional, plane-wave-like 
Bloch orbitals which are completely delocalized over the 
whole crystal £&2d£JLi2&. Wannier functions have been 
used frequently in formal deductions. However, only re- 
cently do they become practically important, e.g., they 
are used to achieve linear scaling with the system size of 
the tight-binding or Kohn-Sham method. 9 ' 10 Moreover, 
as soon as electron correlations in the ground state and in 
excited states are regarded, a local representation offers 
big advantages^i^iiii^^i!^ 

Conventionally, Wannier functions are determined a 
posteriori from Bloch functions utilizing the Wannier 
transformation. 1 ^'^! 6 -' 7 - It is unique up to a unitary ma- 
trix, which can be chosen freely to obtain a certain set of 
Wannier functions, e.g., those which are maximally local- 
ized with respect to a given criterion. Marzari and Van- 
derbilt 6 and Zicovich- Wilson et al£ advocate the Foster- 
Boys functional. 19 ! 20 ' 21 However, there are a number of 
other localization criteria such as the method of Edmis- 
ton and Ruedenberg 2 ^ or the recipe of Pipek and Mezey 23 
which can be used, too. 

Instead, Shukla et al. 2A proposed an embedded clus- 
ter model that facilitates solving modified Hartree-Fock 
equations, which directly yield Wannier-type functions. 
The theory has been implemented in terms of the com- 
puter program WANNIER ) 24 ' 25 which has been used in 
a series of studies of polymers and crystals . 11 ' 24 ' 25 ' 26 ' 27 
It was found, empirically, to provide accurate Hartree- 



Fock energies and band structures. 24 ' 25 ' 27 ' 28 The program 
has also been used as a starting point for post Hartree- 
Fock calculations; the correlation energy of the ground 
and the quasiparticle band structure of several 
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,11,29 



polymers and crystals have been investigated . 12 ' 14 ' 15 ' 17 
The idea to directly solve for Wannier-type functions has 
also been discussed, e.g., in Refs. HHMtand references 
therein) . 

The paper is structured as follows. In Sec. HH I devise 
a rigorous proof of the modified Hartree-Fock equations 
of Shukla et al£^ and give the corresponding modified 
Roothaan-Hall equations. The functions that result from 
these equations are termed a priori Wannier functions. 
The relation of the modified equations to the canonical, 
Bloch-function-based equations and a comparison of the 
computational effort of both approaches can be found 
in Sec. IIIII Conclusions are drawn in Sec. IIV[ specifi- 
cally, I show that the Kohn-Sham equations of density- 
functional theory can be modified in the same way as 
the Hartree-Fock equations to obtain a priori Wannier 
functions. 

Atomic units are used throughout the article. 



II. A PRIORI WANNIER FUNCTIONS 



Definitions 



The sets of functions {w^ (rs), g — 1, . . . , K} are as- 
sociated with each unit cell R of a crystal; they depend 
on the spatial and spin coordinates r and s, respectively. 
The sets are translationally related, i.e., two such sets can 
be brought into coincidence by displacing the functions 
in one of them by a suitable lattice vector R', 

WR+R'eirs) = (rs\R + R'g) =w jig {r-R' s) 

= T R ,w Ag {rs) . (1) 
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Here, Tg, denotes the translation operator for a passive 
translation by R' . Furthermore, the functions are as- 
sumed to be orthonormal with respect to an integration 
over the entire space, 
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The two properties |T]) and @ characterize (spin) Wan- 
nier functionsji^i 3 -^ which are frequently termed (spin) 
Wannier orbitals, if they are Hartree-Fock or Kohn-Sham 
(spin) orbitals. 

I describe the crystal by a nonrelativistic Hamilto- 
nian 5 -i2i 



assuming fixed nuclei. Here, N denotes the number of 
electrons in the crystal which is represented by a paral- 
lelepiped that consists of Nq unit cells. The distance 
between the mth and the nth electron is represented 
by \r m — v n | . The number of nuclei per unit cell is indi- 
cated by M and Z^ A = Za stands for the charge of nu- 



cleus A in unit cell R. Then, 



RA 



is the distance be- 



tween the nth electron and the Ath nucleus in unit cell R. 
Finally, \r^ A — r^, B \ denotes the distance between nu- 
clei A and B of charge Zg A = Za and Zg, B = Zb in 

unit cells R and R' , respectively. 
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B. Wannier-Hartree-Fock equations 



The ansatz for the Hartree-Fock wave function is given 
by the Slater determinant in terms of the occupied Wan- 
nier orbitals, 



N N' 



$o(r 1 s 1 ,...,r N s N ) = A [[[['"s.^fi+li-P' s j+(i-i)JV') i 



(4) 



where N' — is the number of occupied orbitals per unit cell and A stands for the antisymmetrizerj 3 ^ The energy 
expectation value of the Slater determinant ^ with the Hamiltonian reads 



24,32,33 
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it is a functional of the Wannier orbitals. 

In order to minimize the energy expectation value to 
obtain the best description of the ground state by a sin- 
gle Slater determinant in terms of Rltz' variational prin- 
ciple, functional variation with respect to the orbitals in 
Eq. ([5|) is carried out i 14 ' 32 ^ 33 In doing so, I subtract the 
constraints 

N' 

E E A R'SRK(( W R K \ W R'i) ~ S RR' S ^) ( 6 ) 
R,R< 

from E[&q (fi si, . . . , rjv sjv)], employing the Lagrangian 
multipliers A^, ^ this ensures the orthonormality of 
the orbitals. One arrives at the Hartree-Fock equations 



in terms of Wannier orbitals, 

N' 

f\ W R K ) = E A RtR K 1^) 

^ * (7) 

+ E E A R> iR - K » 

with / denoting the Fock operator. Expression ([7]) forms 
a set of N equations which couple the occupied Wannier 
orbitals in a unit cell to the occupied Wannier orbitals 
in all other cells of the crystal. Above all, the Hartree- 
Fock equations ([7]) do not have the form of an eigenvalue 
equation. For these reasons, their practical application 
is cumbersome. 
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To make progress towards a more favorable, modified 
form of Hartree-Fock equations for Wannier orbitals, I 
omit the intercell Lagrangian multipliers in Eq. ([6]) and 
thus arrive at the new energy functional 

C[$$(f*i S!,...,r N s N )] = E[$$(?i S!,...,r N s N )} 

N' 

-EE A ?«(( w r„ I w n i ) ■ 

.(8) 

Here, the translational symmetry of A is exploited, i.e., 
= A(5^jj re = A^ ? ^ K for all lattice vectors R and 
orbital indices £, n. This symmetry, however, will be 
broken again in the next paragraph. Minimizing C leads 
to orthonormal spin-orbitals in each unit cell. However, 
by this simplification of Eq. |6]), one does not enforce 
the mutual intercell orthogonality of the orbitals. Start- 
ing from a properly orthonormalized initial guess for the 
Wannier orbitals in the origin cell may lead to overlap- 
ping orbitals. 

In order to nonetheless achieve orthogonality of the 
occupied orbitals in the origin cell to all other orbitals, 
i.e., their periodic images in all the other unit cells, I 
modify the energy functional C in Eq. (|8|) once more by 
adding an artificial orthogonalizing potential, 

C'[$o(r 1 si,...,r N s N )]= (q) 
C[^(f f 1 s 1 ,...,r N s N )] + Vorth, { ' 

which is defined by 
A N ' 

Vbrth=2 Yl E ( W R K \ W R'i )( W R^\ W Rk) > ( 10 ) 

R^R' 

where A > is called the orthogonalizing potential 
strength or shift parameter. All terms in definition (fT0|) 
are real and non-negative, causing an increase of energy 
proportional to the square of the modulus of the over- 
lap between a pair of orbitals in different unit cells R 
and R' . Minimizing in the limit A — » oo, both pre- 
serves the orthogonality among the Wannier orbitals in 
all cells of the crystal and minimizes the Hartree-Fock en- 
ergy functional E in Eq. ([5]) . It has been shown in practi- 
cal computations that a finite orthogonalization potential 
strength A in the range of 10 3 -I0 5 hartrees causes the re- 
sulting Hartree-Fock energies not to show a noticeable 
dependence on Ai2& 

Expression (JTUJ) can be rewritten compactly by rec- 
ognizing that the orbitals within a particular unit cell 
are orthonormal by construction ([5]). Therefore, the con- 
straint Rk ^ R' £ can be used instead for the summa- 
tion in Eq. (|I0[) . Consequently, the potential Vbrth can 
be expressed in terms of the off-diagonal elements of the 
overlap matrix between the orbitals, 

= (! - S RR> S «t) ( W R K \ W R>t ) • ( n ) 

By taking the trace of if, I arrive at 

Vbrth = ^Tr£ 2 , (12) 



which is an alternate form for Eq. (|10p. 

To minimize the functional £ in Eq. ([9]), I carry out 
functional variation with respect to the N' occupied or- 
bitals in unit cell R^ I arrive at modified Hartree-Fock 
equations 

N' 

(f + ^S)\^R K ) = 12 A RCRn\ W A i ) > (13) 

which I term Wannier- Hartree-Fock equations. The 
penalty projection operator— therein is defined by 

N' 

Vr = E Ek^x^? • ( 14 ) 

R'^rA =1 

It is not translationally symmetric; instead, the rela- 
tion TgVfi = r Pa_g holds. This property of the projec- 
tor Vfi distinguishes the Wannier orbitals in unit cell R 
in the Wannier-Hartree-Fock equations (fl3|) from their 
periodic images in other unit cells because it breaks the 
translational symmetry of the Fock operator /. 

The translational relation of the Wannier function ^ 
implies that it is sufficient to formulate and solve Eq. (fl3|) 
only in the origin cell. For orthogonal orbitals, the La- 
grangian multipliers Ag K constitute a Hermitian ma- 
trix,—^ which is diagonalizable by a unitary transforma- 
tion X^ AX = e. Both / and Vq are invariant under such 
a transformation which mixes the orbitals within every 
unit cell including, particularly, the origin cell. There- 
fore, I formally get a Hermitian N' x N' eigenvalue prob- 
lem 

(/ + AP ff )|%J= EgK \w s J , (15) 

which resembles the canonical Hartree-Fock equa- 
tion o 31 ' 32 i 33 ' 34 i 35 and is thus named pseudocanonical 
Wannier-Hartree-Fock equations. The transformation of 
the orbitals is indicated by affixing a check accent. They 
are referred to as pseudocanonical Wannier orbitals. Due 
to the fact that these orbitals diagonalize the Fock ma- 
trix in the origin cell, they are uniquely determined, apart 
from degeneracies, in analogy to canonical orbitals . 32 ! 33 
Once a self-consistent solution of Eq. (fT5|) has been found, 
the parametrical dependence on A of the orbitals and 
eigenvalues therein vanishes because they are equal to 
the orthonormal orbitals from the previous iteration that 
have been used to construct / and Vq to begin with. 

The pseudocanonical Wannier orbitals are delocalizcd 
over the entire origin cell and thus implicate a similar dis- 
advantageous nonlocality associated with Bloch orbitals 
when applying cutoff criteria to the Fock matrix and to 
the two-electron integrals . 14 ! 15 ' 17 However, this form of 
the Wannier-Hartree-Fock equations is a good starting 
point for further improvements by means of an additional 
localizing potential that can be introduced in the expres- 
sion for C' in Eq. (J9j> . Established forms of localizing 
potentials are the one of Edmiston and Ruedenberg^ 2 - for 
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their localization criterion and the one of Gilbert^ for 
the Foster-Boys criterion! 19 i 20 i 21 

The functional dependence of / + A Va on the occupied 
Wannier orbitals can be disregarded) 32 ' Then, the Fock 
operator becomes a conventional Hermitian operator and 
the restriction of the penalty projection operator, to act 
only on occupied orbitals, can be released, i.e., N' in 
Eq. (fT4|) is replaced by the total number of orbitals per 
unit cell K. Now, Eq. (fTS)) holds also for virtual Wannier 
orbitals. 



C. Wannier-Roothaan-Hall equations 

Spin Wannier orbitals have been employed so far. Let 
me assume a restricted, closed-shell Hartree-Fock point 
of view to remove the spin dependence! 32 ' 33 In this case, 
the spin orbitals are expressed as the product of a spatial 
orbital w^ g (r) with the spinor for spin up a(s) and spin 
down (3(s), respectively. The spatial orbitals, wj^ g (r), 
q = 1, . . . ,K, are expanded in terms of one-particle basis 
functions,- 7 -^ Xn{r), M = 1, ■ • • ,K, 

k 



g M=l 



S 2 ^g+R u RqXh{t 



9-R) 



(16) 

The number of basis functions K determines the num- 
ber of orbitals per unit cell. In Eq. (fT6|) . I exploit the 
fact that both the spatial Wannier orbitals and the ba- 
sis functions Tg + g XmM f° rm se ts of functions whose 
members are related by lattice translations ([1]). The ex- 
pansion coefficients, hence, are translationally symmet- 
ric C g+Rn$, e = C 3n o e = C ne(9)- The displacement d M 
of the /xth basis function in a unit cell accounts for the 
fact that a basis function is frequently centered on atoms 
which, in turn, are displaced somewhat from the origin 
of the unit cell. 

The basis-set representation of the Fock opera- 
tor / is Fgpfv = (Xg^\f\Xg'u), the overlap ma- 
trix is Sg-^g-'y = (Xgfj, I Xg'v)-, and the matrix of 
the penalty projection operator (|14p — with TV' replaced 
by K — reads 2 ^ 



^Rlgng'v ~ ( Xgy, I ^r I Xg" u ) 



R',91,92 K i?>e=i 

x C g t kR'b C ~ 2 ( r, B S 3' vgii 



(17) 



Representing Eq. (fT5| in the basis set, I obtain modified 
Roothaan-Hall equations 2 ^ 



(F + XV 5 )C B = e 5 SC e 



(18) 



to which I refer as Wannier-Roothaan-Hall equations. 
Formula (|18|) yields the occupied and virtual Wannier 



orbitals in the origin cell (r) with g = 1, . . . , K and 
expansion coefficients (Cg)^^ = C g . fi tf g for all K Nq in- 
dices gfi. In expression (fT8|) . the coupling of the Wan- 
nier orbitals in the origin cell to the Wannier orbitals in 
neighboring unit cells reappears by means of the basis-set 
representation ()16j) which could be avoided by omitting 
the Lagrangian multipliers in Eq. (jSJ). Consequently, the 
dimension of the matrices in Eq. (|18[) scales both with 
the number of basis functions in the origin cell and with 
the number of unit cells utilized to support the Wannier 
orbitals. Note that only a subset of K eigenvectors out of 
the K No eigenvectors of Eq. (TT5)) is required. The Wan- 
nier orbitals in other but the origin cell are simply given 
by lattice translations (TIJ, exploiting that the expansion 
coefficients C^ g (g) are translationally symmetric. There- 
with, all Wannier orbitals of the crystal are determined. 
The energy bands are not given by the e - g in Eq. (| 18|) . 
Instead, they are found by diagonalizing the Hermitian 

matrix F(k) = Yl elkR -^or f° r a number of k points, 
R 

where (F^^) ae = (w^ a \ f ) are blocks of the Fock 
matrix in Wannier representation. 

The solution of the Wannier-Roothaan-Hall equa- 
tions (fT8|) has been implemented in the wannier pro- 
gra m 24 i 25 and its accuracy was tested in a series of 
studies. Systematic comparisons were made of ground- 
state properties and band structures of crystal o 25 ' 27 
and infinite chains^ The basis-set expansion (IT51) is 
found to provide a satisfactory accuracy using up to 
third-nearest-neighbor cells for ionic three-dimensional 
crystal o 14 ' 15 ' 24 ' 25 ' 26 ' 27 as well as covalently bonded and 
hydrogen-bonded infinite chains . 13 i 14 i 16 i 17 i 28 Note, how- 
ever, that despite the fact that variation of the expansion 
coefficients in Eq. (fT6|) is allowed only in a small clus- 
ter of unit cells, a sophisticated treatment of the very 
long range electrostatic interactions in the Fock operator 
is required, accounting for a large number of unit cells. 
Therefore, one uses the Ewald summation technique for 
the Coulomb matrix elements and an extended summa- 
tion for the exchange matrix elements^ The occupied 
pseudocanonical a priori Wannier orbitals are found to 
reach a few angstroms from the atom they are centered 
on» 



III. PROPERTIES 



A. Relation to canonical Hartree-Fock equations 

The connection of the Wannier-Hartree-Fock equa- 
tions (TT5]) to the corresponding equations of a Bloch- 
orbital-based formalism — the canonical Hartree-Fock 
equation o 32 ' 33 — can be achieved readily. To this end, 
the orthogonalizing potential (112|) is replaced by an ex- 
pression which comprises the full overlap matrix, i.e., 
^ Tr S 2 . As the Wannier functions are normalized to 

unity, one obtains Tr S = Tr S, + K Nq. This amounts 
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to a meaningless overall energetic shift of £' in Eq. ([9]) 
upon replacing Vo r th- Functional variation of the new £ 
yields that the summation in the penalty projector of 
Eqs. I|14p and (|17[) is no longer restricted to the unit 
cells in the neighborhood but also contains the origin 
cell, i.e., the operator 1S exchanged by The 

R'^R R' 

new projector V becomes translationally symmetric and 
actually the identity operator 1. This transformation 
of the Hartree-Fock equations is equivalent to the addi- 

N' 

tion of I w Re ) ( w Re I on both sides of Eq. (fT3|) . which 

causes the eigenvalues e$ K [cf. Eq. (fT5|) ] to shift by A, 
i.e., they become e^ K + A. The modification undoes the 
orthogonalizing potential and reverts the equations to 
the modified Hartree-Fock equations which result from 
Eq. ((8]) with its eigenvalues shifted by A. Expression ([7|) 
is obtained again by realizing that the matrix representa- 
tion of / + A V in terms of Wannier orbitals also contains 
off-diagonal terms with respect to the lattice vectors, thus 
reintroducing off-diagonal Lagrangian multipliers. 

The arguments of the previous paragraph can be ex- 
pressed more clearly by changing the Roothaan-Hall 
equations (|18|) because their dimensionality and formal 
structure are preserved under the replacement of the pro- 
jector. Due to the translational symmetry of / + V, 
Born-von Karman boundary conditions become benefi- 
ciali^ With the matrix representation of V = I, which 
is V = S, the modified Roothaan-Hall equations (fT8|) are 
changed tc-^ 

(F + XV) C = SC'e' . (19) 

The N K x N K matrices F, S, and V are cyclic matri- 
ces which can be block diagonalized employing the uni- 
tary transformational W^ g ^ p = -^=8 ep e lkR . Mul- 
tiplying with from the left and inserting VVW^ be- 
fore C , the eigenvalue problem (p~9j) breaks down into 
No independent K x K eigenvalue problems, 

[F(k) + XV(k)] C'(k) = S(k)C(k)e'(k) . (20) 

Solving these equations yields orthonormal Bloch orbitals 
and energy bands which are shifted by A. 



B. Computational efficiency 

To compare the computational effort of the canoni- 
cal Hartree-Fock method to the modified Hartree-Fock 
method, I analyze the corresponding Roothaan-Hall 
equations, Eqs. (|2"01) and (fT5)l . respectively. Without re- 
garding space-group symmetry, both sets of equations 
have the same dimensionality because the number of or- 
bitals (or basis functions) is the same, namely, K Nq. 
The effort to compute the matrix elements is clearly the 
same for the basis-set overlap matrix; it is also equal 



for the Fock matrix because the Fock operator is invari- 
ant under unitary transformations of the orbitals. The 
support of the basis-set expansion is the same for Bloch 
and Wannier orbitals: the basis functions in the paral- 
lelepiped of No unit cells. Only the expansion coefficients 
differ. Consequently, the number of floating point oper- 
ations necessary to determine the Fock operator is the 
same, if no further approximations are made. Yet the 
projection operator, which is only required for the mod- 
ified Hartree-Fock equations, requires an extra effort. 

In a next step, the modified and canonical Roothaan- 
Hall equations need to be diagonalized. In both cases, 
a full diagonalization requires the same effort. How- 
ever, in practice a selective computation is carried out. 
The canonical equations are block diagonalized first, 
Eq. (j2"0|) , and the subblocks arc independently diago- 
nalized fully afterwards for a grid of k points^ The 
Wannier- Roothaan-Hall equations have to be treated dif- 
ferently. The spectrum of F + APg in Eq. (fTB)) has a very 
favorable property: the lower K eigenvalues correspond 
to the Wannier orbitals in the origin cell. The other 
eigenvalues are well separated from the former ones be- 
cause they are shifted to high values by the shift parame- 
ter Ai2£ Iterative eigenvalue solvers, particularly the one 
of Davidson^ can be employed to reduce the numerical 
effort to determine the lowest K eigenpairs of F + A'Pg. 
Their effort is predominantly determined by the matrix 
times vector product. As the summation expression for 
this product and the summations for the block diago- 
nalization and subsequent fc-space integration are sim- 
ilar, the overall computational effort of both methods, 
to determine crystal orbitals with a certain accuracy, 
should be comparable. In fact, the Bloch-orbital-based 
equations can be solved more efficiently for crystals with 
a small number of atoms per unit cell; the solution of 
the Wannier-orbital-based equations is more efficient for 
crystals with a large unit cell because cutoff criteria can 
be established to lower the actual effort. 



IV. CONCLUSION 

In this paper, I derive modified Hartree-Fock equa- 
tions, which directly yield Wannier orbitals. An orthog- 
onalizing potential is added to the Hartree-Fock energy 
functional to ensure the proper orthonormality of the 
resulting orbitals. It serves to replace the Lagrangian 
multipliers between unit cells needed otherwise. The 
equations necessarily break the translational symmetry 
as Wannier functions are translationally related, in con- 
trast to Bloch functions which are translationally sym- 
metric. I show how the conventional Bloch-orbital-based 
Hartree-Fock equations can be recovered by restoring the 
translational symmetry of the modified equations. Ana- 
lyzing the spectral properties of the modified Fock ma- 
trix, I find the numerical effort of the method to be com- 
parable to the effort of the Bloch-orbital-based approach. 

The orbitals which result from the modified Hartree- 
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Fock equations are referred to as pseudocanonical Wan- 
nier functions as they are delocalized over the entire 
unit cell. However, they can be localized additionally 
within unit cells by adding a suitable localizing potential 
to the energy expression. Particularly, the potential of 
Gilbert 36 is to be mentioned here, which minimizes the 
Foster-Boys functionah 19 ' 20 ' 21 It thus offers a different 
route to Refs.l6||7|to determine maximally localized Wan- 
nier functions. Alternatively, the potential of Edmiston 
and Rucdcnberg can be used as a localizing potential* 2 - 2 * 
This study is based on the Hartree-Fock theory. How- 
ever, the extension of the ideas to density-functional 
theor y 40 i 41 ' 42 is straightforward. To this end, one ex- 
changes the Hartree-Fock energy functional (O by the 
Hohenberg-Kohn energy functional in terms of the Kohn- 
Sham orbitalsi 41 i 42 The Hohenberg-Kohn variational the- 
ore m 40 ' 42 ensures that minimizing this functional with 
respect to the orbitals provides the exact ground-state 
energy, if one uses the exact exchange-correlation energy 
functional. The minimization is carried out under the 



constraint that the orbitals remain orthonormal using 
Eq. ([6|). In other words, we can essentially follow the line 
of argument that leads from Eq. ([5]) to Eq. (fT"5|) . I term 
the latter equation — with the Fock operator replaced 
by the Kohn-Sham operator — pseudocanonical Wannier- 
Kohn-Sham equations. All subsequent arguments and 
expressions carry over analogously. 
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